import math
import typing as tp
from math import (
    pi as PI,
    sin, atan,
    sqrt,
    fabs, cos,
)

import numpy as np
from matplotlib import pyplot as plt

from utils import ln, degree2radian

A_ELEV = PI / 8
BW = 15 / 180 * PI
R_M = 1000


# 干扰机带宽 jammer bandwidth;
# 干扰机功率 jammer power;
# 干扰机距离 jammer range;
# 干扰机探测雷达 jammer finder;
# 干扰机天线增益 jammer antenna gain;
# 干扰机引导设备 jammer directing equipment

class RadarParameter(object):
    def __init__(self):
        # P_t 为雷达天线的发射功率;
        # P_r 为雷达天线的接收功率; 通常与发射功率相等（单基地雷达）

        # G_t 为雷达发射天线功率增益;
        # G_r 为雷达接收天线功率增益;

        # sigma 为雷达散射截面积RCS;
        # lambda_ 为雷达发射波长;

        # k 为玻尔兹曼常数;
        # F_n 是噪声系数;
        # min_S_N 为最小检测信噪比;
        # T_0 为接收机噪声温度;
        # B_n 为接收机同频带宽;
        # L 是系统损耗因子。
        self.P_t = 10 ** 6  # 发射功率/W
        self.P_r = 10 ** 6
        self.G_t = 40  # 天线增益/dB
        self.G_r = 40

        self.sigma = 3  # 目标截面积/m2
        self.lambda_ = 0.056  # 波长/m

        self.n = 16  # 脉冲积累数

        self.k = 1.38e-23
        self.F_n = 1
        self.min_S_N = 2  # 最小检测信噪比
        self.T_0 = 1
        self.B_n = 1600000  # 接收机带宽/Hz
        self.L = 10  # 系统损耗因子/dB

        self.vertical_bw = degree2radian(5)  # 雷达垂直波束宽度/°
        self.horizontal_bw = degree2radian(5)  # 雷达水平波束宽度/°

        self.A_ELEVATION = degree2radian(30)  # 天线仰角/°(雷达仰角？)


class JammerParameter(object):
    K = 2  # 压制系数

    def __init__(self):
        pass
        # K 为干扰压制系数;
        # n 为脉冲积累系数;

        # Pi 为干扰机 i 的干扰功率;
        # Bi 为干扰机 i 的发射带宽;
        # γi为干扰机 i 的干扰信号的极化损耗;

        # Li 为干扰机 i 的系统损耗因子;
        # Gi为干扰机 i 在雷达方向上的增益。

        # self.K = 2  # 压制系数

        self.P_i = 1  # 干扰功率/W
        self.G_i = 30  # 天线增益/dB
        self.gamma_i = 0.5  # 极化损耗

        # 主瓣宽度/° 10

        self.R_i = 100  # Ｒi 为干扰机 i 到雷达天线的距离;
        self.B_i = 2000000  # 发射带宽/Hz
        self.L_i = 10  # 系统损耗因子/dB

        self.xi_i = degree2radian(300)  # ξi干扰机在局部坐标系中方位角为
        self.delta_i = degree2radian(60)  # 天线俯仰角/°
        self.K_gain = 0.08  # 增益修正系数


def a_f(bw):
    return 1.39157 / sin(bw / 2)


A_F = a_f(BW)


def a_max(a_f):
    # -PI/2, PI/2
    return atan(PI / (
            a_f * sqrt(1 - (PI / a_f) ** 2)
    ))


def F_theta_phi_without_jammer_sinx(theta, phi, radar_parameter: RadarParameter):
    # radar antenna direction
    """雷达天线方向图：(sin x)/x 型"""
    delta_a = phi - radar_parameter.A_ELEVATION
    if fabs(delta_a) <= 1e-6:
        return 1
    if fabs(delta_a) > a_max(A_F):
        return 0.03
    t = A_F * sin(delta_a)
    return sin(t) / t


def F_theta_phi_without_jammer_gaussian2(theta, phi, radar_parameter: RadarParameter):
    theta_b = radar_parameter.vertical_bw
    k = 4 * ln(math.sqrt(2)) * (theta_b) ** 2
    f = math.e ** (-k * (phi ** 2))
    return f


def F_theta_phi_without_jammer_gaussian(theta, phi, radar_parameter: RadarParameter):
    """
    theta
    phi

    theta_b: 雷达波束宽度，
    phi_z: 雷达仰角

    :return:
    :rtype:
    """
    theta_b = radar_parameter.vertical_bw
    phi_z = radar_parameter.A_ELEVATION
    k_0 = ln(2) / (2 * (math.sin(theta_b) ** 2))
    _ = math.sin(phi) - math.sin(phi_z)
    # return math.e ** (-2.766 * (_ ** 2)) / ((phi_z / 2) ** 2)
    return math.e ** (k_0 * (_ ** 2))


# 通用的公式
def _coverage(theta, phi, R_max, F):
    return R_max(theta, phi) * F(theta, phi)


def coverage_without_jammer_sinx(theta, phi):
    return R_M * F_theta_phi_without_jammer_sinx(theta)


def __coverage_ideal(theta, phi):
    """
    theta 为方位角; phi 为俯仰角;

    sigma 为雷达散射截面积;
    P_t 为雷达天线的发射功率;
    G_r 为雷达接收天线功率增益;
    F_r 为雷达接收端到目标的方向图因子;
    lambda_ 为雷达发射波长;
    G_t 为雷达发射天线功率增益;
    F_t 为雷达发射端到目标的方向图因子;
    k 为玻尔兹曼常数;
    F_n 是噪声系数;
    min_S_N 为最小检测信噪比;
    T_0 为接收机噪声温度;
    B_n 为接收机同频带宽;
    L 是系统损耗因子。
    :param theta:
    :type theta:
    :param phi:
    :type phi:
    :return:
    :rtype:
    """
    sigma = 3
    P_t = 10 ** 6
    G_r = 40
    F_r = F_theta_phi_without_jammer_gaussian(theta, phi)
    lambda_ = 0.056
    G_t = 40
    F_t = F_theta_phi_without_jammer_gaussian(theta, phi)
    k = 1.38e-23
    F_n = 1
    min_S_N = 2
    T_0 = 1
    B_n = 1.6 * 10e6
    L = 10

    _1 = (
            P_t * G_r * G_t * sigma *
            (lambda_ * F_r * F_t) ** 2
    )
    _2 = (
            ((4 * math.pi) ** 3)
            * k
            * T_0
            * B_n
            * F_n
            * min_S_N
            * L
    )
    _ = _1 / _2

    return _ ** 0.25


def jammer_value(theta, phi, radar_parameter: RadarParameter, jammer_parameter: JammerParameter):
    # theta = degree2radian(10)
    if 0 <= theta <= 0.5 * radar_parameter.horizontal_bw:
        return radar_parameter.G_r
    else:
        # 0.5 * radar_parameter.theta_bw < theta < 0.5 * math.pi:
        return jammer_parameter.K_gain * (
                radar_parameter.horizontal_bw / min(theta, 0.5 * math.pi)) * radar_parameter.G_r


def R_max_with_jammer(theta, phi, radar_parameter: RadarParameter, jammer_parameter: JammerParameter):
    _1 = (
            1
            * radar_parameter.P_t
            * (radar_parameter.G_t ** 2)
            * jammer_parameter.K
            * radar_parameter.sigma
            * (radar_parameter.lambda_ ** 2)
            * (radar_parameter.n ** 0.5)
    )
    abs_theta_xi = abs(theta - jammer_parameter.xi_i)
    if abs_theta_xi <= math.pi:
        jammer_value_theta = abs_theta_xi
    else:
        jammer_value_theta = 2 * math.pi - abs_theta_xi

    _2 = (
            4 * math.pi
            * jammer_parameter.P_i
            * jammer_parameter.G_i
            * jammer_value(jammer_value_theta, None, radar_parameter, jammer_parameter)
            * radar_parameter.B_n
            * jammer_parameter.gamma_i
            / (
                    (jammer_parameter.R_i ** 2) * jammer_parameter.B_i * jammer_parameter.L_i
            )
    )

    return (_1 / _2) ** 0.25


def R_max_with_jammers(theta, phi, radar_parameter: RadarParameter, jammer_parameters: tp.List[JammerParameter]):
    assert len(jammer_parameters) > 0

    _1 = (
            1
            * radar_parameter.P_t
            * (radar_parameter.G_t ** 2)
            * jammer_parameters[0].K
            * radar_parameter.sigma
            * (radar_parameter.lambda_ ** 2)
            * (radar_parameter.n ** 0.5)
    )

    _2_sum = 0
    for jammer_parameter in jammer_parameters:
        abs_theta_xi = abs(theta - jammer_parameter.xi_i)
        if abs_theta_xi <= math.pi:
            jammer_value_theta = abs_theta_xi
        else:
            jammer_value_theta = 2 * math.pi - abs_theta_xi

        _2_sum += (1
                   * jammer_parameter.P_i
                   * jammer_parameter.G_i
                   * jammer_value(jammer_value_theta, None, radar_parameter, jammer_parameter)
                   * radar_parameter.B_n
                   * jammer_parameter.gamma_i
                   / (
                           (jammer_parameter.R_i ** 2) * jammer_parameter.B_i * jammer_parameter.L_i
                   )
                   )
    _2 = 4 * math.pi * _2_sum

    return (_1 / _2) ** 0.25


def R_max_without_jammer(theta, phi):
    return R_M


# receive, transmit
def F_theta_phi_with_jammer_gaussian(theta, phi, radar_parameter: RadarParameter, jammer_parameter: JammerParameter):
    abs_delta_phi = abs(jammer_parameter.delta_i - phi)
    if abs_delta_phi <= 0.5 * radar_parameter.vertical_bw:
        return F_theta_phi_without_jammer_gaussian(theta, phi, radar_parameter)
    else:
        return (
                1 / jammer_parameter.K_gain
                * (jammer_parameter.delta_i / radar_parameter.vertical_bw)
                * F_theta_phi_without_jammer_gaussian(theta, phi, radar_parameter)
        )


# receive, transmit
def F_theta_phi_with_jammer_gaussian2(theta, phi, radar_parameter: RadarParameter, jammer_parameter: JammerParameter):
    abs_delta_phi = abs(jammer_parameter.delta_i - phi)
    if abs_delta_phi <= 0.5 * radar_parameter.vertical_bw:
        return F_theta_phi_without_jammer_gaussian2(theta, phi, radar_parameter)
    else:
        return (
                1 / jammer_parameter.K_gain
                * (jammer_parameter.delta_i / radar_parameter.vertical_bw)
                * F_theta_phi_without_jammer_gaussian2(theta, phi, radar_parameter)
        )


def F_theta_phi_with_jammer_sinx(theta, phi, radar_parameter: RadarParameter, jammer_parameter: JammerParameter):
    abs_delta_phi = abs(jammer_parameter.delta_i - phi)
    if abs_delta_phi <= 0.5 * radar_parameter.vertical_bw:
        return F_theta_phi_without_jammer_sinx(theta, phi, radar_parameter)
    else:
        return (
                1 / jammer_parameter.K_gain
                * (jammer_parameter.delta_i / radar_parameter.vertical_bw)
                * F_theta_phi_without_jammer_sinx(theta, phi, radar_parameter)
        )


def coverage_with_jammer_gaussian(theta, phi, radar_parameter: RadarParameter, jammer_parameter: JammerParameter):
    return (
            R_max_with_jammer(theta, phi, radar_parameter, jammer_parameter)
            * (F_theta_phi_with_jammer_gaussian(theta, phi, radar_parameter, jammer_parameter) ** 0.5)
    )


def coverage_with_jammer_gaussian2(theta, phi, radar_parameter: RadarParameter, jammer_parameter: JammerParameter):
    return (
            R_max_with_jammer(theta, phi, radar_parameter, jammer_parameter)
            * (F_theta_phi_with_jammer_gaussian2(theta, phi, radar_parameter, jammer_parameter) ** 0.5)
    )


def coverage_with_jammers_gaussian2(theta, phi, radar_parameter: RadarParameter,
                                    jammer_parameters: tp.List[JammerParameter]):
    assert len(jammer_parameters) > 0
    return (
            R_max_with_jammers(theta, phi, radar_parameter, jammer_parameters)
            * (F_theta_phi_with_jammer_gaussian2(theta, phi, radar_parameter, jammer_parameters[0]) ** 0.5)
    )


def coverage_with_jammer_sinx(theta, phi, radar_parameter: RadarParameter, jammer_parameter: JammerParameter):
    return (
            R_max_with_jammer(theta, phi, radar_parameter, jammer_parameter)
            * (F_theta_phi_with_jammer_sinx(theta, phi, radar_parameter, jammer_parameter) ** 0.5)
    )


# def coverage_with_jammer_sinx(theta, phi, radar_parameter: RadarParameter, jammer_parameter: JammerParameter):
#     return (
#             R_max_with_jammer(theta, phi, radar_parameter, jammer_parameter)
#             * F_theta_phi_with_jammer_sinx(theta, phi, radar_parameter, jammer_parameter)
#     )


def sampling(coverage):
    x, y, z = [], [], []

    step_theta, step_phi = 0.1, 0.02
    theta_min, theta_max = [math.pi * _ for _ in [0, 2]]
    theta = theta_min
    # ，θ为雷达天线的方位角，φ为雷达天线的俯仰角
    while theta <= theta_max:
        _sin_theta = sin(theta)
        _cos_theta = cos(theta)

        phi_min, phi_max = [math.pi * _ for _ in [-0.5, 0.5]]
        # phi_min, phi_max = [math.pi * _ for _ in [0, 1]]

        phi = phi_min
        while phi <= phi_max:
            _sin_phi = sin(phi)
            _cos_phi = cos(phi)

            R = coverage(theta, phi)
            x.append(R * _cos_theta * _cos_phi)
            y.append(R * _sin_theta * _cos_phi)
            z.append(R * _sin_phi)

            phi += step_phi

        theta += step_theta

    print(len(x))
    return x, y, z


def plot_coverage(x, y, z):
    x, y, z = np.array(x), np.array(y), np.array(z)
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    # x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
    # y = [5, 6, 7, 8, 2, 5, 6, 3, 7, 2]
    # z = [1, 2, 6, 3, 2, 7, 3, 3, 7, 2]
    # X, Y = np.meshgrid(x, y)
    # Z = X * np.exp(X - Y)
    # ax.plot_surface(X, Y, Z)

    ax.scatter3D(x, y, z, cmap='Blues', s=2)  # 绘制散点图
    ax.plot3D(x, y, z, 'gray', linewidth=0.5, markersize=12)  # 绘制空间曲线

    # fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
    # surf = ax.plot_surface(x, y, z, cmap=cm.coolwarm,
    #                        linewidth=0, antialiased=False)

    plt.show()


# def save(points)
#     pass

if __name__ == '__main__':
    coverage = coverage_without_jammer_sinx

    radar_parameter = RadarParameter()
    jammer_parameter = JammerParameter()
    # for i in range(10):
    #     xi = 36 * i
    #     print(xi)
    #     jammer_parameter.xi_i = degree2radian(xi)
    #     coverage = lambda theta, phi: coverage_with_jammer_gaussian2(theta, phi, radar_parameter, jammer_parameter)
    #     r = sampling(coverage)
    #     plot_coverage(*r)

    # coverage = lambda theta, phi: coverage_with_jammer_sinx(theta, phi, radar_parameter, jammer_parameter)
    # coverage = lambda theta, phi: coverage_with_jammer_gaussian(theta, phi, radar_parameter, jammer_parameter)
    # coverage = lambda theta, phi: coverage_with_jammer_gaussian2(theta, phi, radar_parameter, jammer_parameter)

    # 多个干扰机
    j1 = JammerParameter(); j1.xi_i = degree2radian(30)
    j2 = JammerParameter(); j2.xi_i = degree2radian(60)
    j3 = JammerParameter(); j3.xi_i = degree2radian(180)
    j4 = JammerParameter(); j3.xi_i = degree2radian(120)
    jammer_parameters = [j1, j2, j3, j4]
    coverage = lambda theta, phi: coverage_with_jammers_gaussian2(theta, phi, radar_parameter, jammer_parameters)

    r = sampling(coverage)
    # # print(r)
    plot_coverage(*r)
